
cap log close	
local gruppo = 1
	
log  using  .../AA_8_`gruppo'_cap_3_5Sep2020.smcl,replace

version 13	
use .../sample_estimation_`gruppo'_rev1_cap_3_Sep2020.dta, clear
		
*** 1 ***	
	cap program drop prog2_lnf
	set type double					

	program define prog2_lnf

	args todo b lnf

	tempvar ISE IU ISELF IULF logitp2 logitp3 logitp4 b1 b2 b3 b4   a1 a2 a3 a4 logitp5 logitp6 logitp7 logitp8  lnfie last lE1 lE2 lE3   lE1LF lE2LF lE2LF ///
	        lU1  lU2 lU3 lU1LF lU2LF lU3LF l1 ll1 l2 ll2 l3 ll3 l4 ll4 l5 ll5 l6 ll6 l7 ll7 l8 ll8  mlSE1 mlSE2  mlSE3  mlE1 mlE2 mlE3 mlSELF1  mlSELF2 mlSELF3  mlELF1 mlELF2 mlELF3 ///
	        h1SE h2SE h3SE  sum1SE sum2SE  sum3SE  ST1SE ST2SE ST3SE lE3LF h1SELF h2SELF  h3SELF  sum1SELF sum2SELF  sum3SELF  ST1SELF ST2SELF ST3SELF ///
		h1U h2U  h3U  sum1U sum2U  sum3U  ST1U ST2U ST3U h1ULF h2ULF   h3ULF  sum1ULF sum2ULF  sum3ULF  ST1ULF ST2ULF ST3ULF lastSE lastU lastSELF lastULF lnfi
		
	tempname p2 p3 p4 p1
	
        *defining equations
    mleval `ISE'     = `b', eq(1)
	mleval `IU'      = `b', eq(2)
    mleval `ISELF'   = `b', eq(3)
	mleval `IULF'    = `b', eq(4)
	mleval `a1'       = `b', scalar eq(5)
	mleval `a2'       = `b', scalar eq(6)
	mleval `b1'       = `b', scalar eq(7)
	mleval `b2'       = `b', scalar eq(8)
	mleval `logitp2' = `b', scalar  eq(9) //5



	

	quietly {
			* to make sure probabilities constrained between 0 and 1
			* for identification, p11==1-p21-p22 and m11==0
			scalar `p2' = exp(`logitp2')/(1+exp(`logitp2'))
			scalar `p1' = 1             /(1+exp(`logitp2'))

			***** for SELFEMPLOYED spells
			
			* hazard functions for the three types
			by $S_E_id $spell: gen double `h1SE' = 1-exp(-exp(`ISE' + `a1'))	if $SE_E ==1 & $LF == 0
			by $S_E_id $spell: gen double `h2SE' = 1-exp(-exp(`ISE' ))			if $SE_E ==1 & $LF == 0
			
                        * first part of survival function 
			by $S_E_id $spell: gen double `sum1SE' = sum(exp(`ISE' + `a1'))		if $SE_E ==1 & $LF == 0
			by $S_E_id $spell: gen double `sum2SE' = sum(exp(`ISE' ))			if $SE_E ==1 & $LF == 0
		
			
			* cumulative product of survival functions
			by $S_E_id $spell: gen double `ST1SE' = exp(-`sum1SE'[_N])			if _n==_N & $SE_E ==1 & $LF == 0
			by $S_E_id $spell: gen double `ST2SE' = exp(-`sum2SE'[_N])			if _n==_N & $SE_E ==1 & $LF == 0
			

			* identify last obs of the spell
			by $S_E_id $spell: gen byte `lastSE' = (_n==_N) if $SE_E ==1 & $LF == 0
			replace `lastSE'=0 if `lastSE'==.

			* likelihood functions of the spells for each type
			by $S_E_id $spell: gen double `lE1' =  `ST1SE' * ((`h1SE'/(1-`h1SE'))^$event ) if _n==_N & $SE_E ==1 & $LF == 0
			by $S_E_id $spell: gen double `lE2' =  `ST2SE' * ((`h2SE'/(1-`h2SE'))^$event ) if _n==_N & $SE_E ==1 & $LF == 0
			
			
			replace `lE1'=0 if `lE1'==.
			replace `lE2'=0 if `lE2'==.
			

			******** for EMPLOYED 
			by $S_E_id $spell: gen double `h1U' = 1-exp(-exp(`IU' + `b1' )) if $SE_E ==0 & $LF == 0
			by $S_E_id $spell: gen double `h2U' = 1-exp(-exp(`IU' )) if $SE_E ==0 & $LF == 0
			
			
			by $S_E_id $spell: gen double `sum1U' = sum(exp(`IU' + `b1' )) if $SE_E ==0 & $LF == 0
			by $S_E_id $spell: gen double `sum2U' = sum(exp(`IU' )) if $SE_E ==0 & $LF == 0
			

			by $S_E_id $spell: gen double `ST1U' = exp(-`sum1U'[_N]) if _n==_N & $SE_E ==0 & $LF == 0
			by $S_E_id $spell: gen double `ST2U' = exp(-`sum2U'[_N]) if _n==_N & $SE_E ==0 & $LF == 0
			
			
			by $S_E_id $spell: gen byte `lastU' = (_n==_N) if $SE_E ==0 & $LF == 0
			replace `lastU'=0 if `lastU'==.

			by $S_E_id $spell: gen double `lU1' =  `ST1U' * ((`h1U'/(1-`h1U'))^$event ) if _n==_N & $SE_E ==0 & $LF == 0
			by $S_E_id $spell: gen double `lU2' =  `ST2U' * ((`h2U'/(1-`h2U'))^$event ) if _n==_N & $SE_E ==0 & $LF == 0
                        
			               
			replace `lU1'=0 if `lU1'==.
			replace `lU2'=0 if `lU2'==.
			
			
			***** for SELFEMPLOYED left-censored spells
			
			* hazard functions for the three types
			by $S_E_id $spell: gen double `h1SELF' = 1-exp(-exp(`ISELF' + `a2'))   if $SE_E ==1 & $LF == 1
			by $S_E_id $spell: gen double `h2SELF' = 1-exp(-exp(`ISELF'))   if $SE_E ==1 & $LF == 1
			
			
                        * first part of survival function
			by $S_E_id $spell: gen double `sum1SELF' = sum(exp(`ISELF' + `a2'))  if $SE_E ==1 & $LF == 1
			by $S_E_id $spell: gen double `sum2SELF' = sum(exp(`ISELF' ))    if $SE_E ==1 & $LF == 1
			
			
			* cumulative product of survival functions
			by $S_E_id $spell: gen double `ST1SELF' = exp(-`sum1SELF'[_N]) if _n==_N & $SE_E ==1 & $LF == 1
			by $S_E_id $spell: gen double `ST2SELF' = exp(-`sum2SELF'[_N]) if _n==_N & $SE_E ==1 & $LF == 1
			
			
			* identify last obs of the spell
			by $S_E_id $spell: gen byte `lastSELF' = (_n==_N) if $SE_E ==1 & $LF == 1
			replace `lastSELF'=0 if `lastSELF'==.

			* likelihood functions of the spells for each type
			by $S_E_id $spell: gen double `lE1LF' =  `ST1SELF' * ((`h1SELF'/(1-`h1SELF'))^$event ) if _n==_N & $SE_E ==1 & $LF == 1
			by $S_E_id $spell: gen double `lE2LF' =  `ST2SELF' * ((`h2SELF'/(1-`h2SELF'))^$event ) if _n==_N & $SE_E ==1 & $LF == 1
			
           
			replace `lE1LF'=0 if `lE1LF'==.
			replace `lE2LF'=0 if `lE2LF'==.
			
			

			******** for EMPLOYED left-censored spells
			by $S_E_id $spell: gen double `h1ULF' = 1-exp(-exp(`IULF' + `b2')) if $SE_E ==0 & $LF == 1
			by $S_E_id $spell: gen double `h2ULF' = 1-exp(-exp(`IULF' )) if $SE_E ==0 & $LF == 1
			

			by $S_E_id $spell: gen double `sum1ULF' = sum(exp(`IULF' + `b2')) if $SE_E ==0 & $LF == 1
			by $S_E_id $spell: gen double `sum2ULF' = sum(exp(`IULF' )) if $SE_E ==0 & $LF == 1
			
			
			by $S_E_id $spell: gen double `ST1ULF' = exp(-`sum1ULF'[_N]) if _n==_N & $SE_E ==0 & $LF == 1
			by $S_E_id $spell: gen double `ST2ULF' = exp(-`sum2ULF'[_N]) if _n==_N & $SE_E ==0 & $LF == 1
			
			
			by $S_E_id $spell: gen byte `lastULF' = (_n==_N) if $SE_E ==0 & $LF == 1
			replace `lastULF'=0 if `lastULF'==.

			by $S_E_id $spell: gen double `lU1LF' =  `ST1ULF' * ((`h1ULF'/(1-`h1ULF'))^$event ) if _n==_N & $SE_E ==0 & $LF == 1
			by $S_E_id $spell: gen double `lU2LF' =  `ST2ULF' * ((`h2ULF'/(1-`h2ULF'))^$event ) if _n==_N & $SE_E ==0 & $LF == 1
			
                                       
			replace `lU1LF'=0 if `lU1LF'==.
			replace `lU2LF'=0 if `lU2LF'==.
			
	    
				// individual log likelyhood for each type
			by $S_E_id: gen double `l1'=`lE1' +`lU1'+`lU1LF'+`lE1LF'
			by $S_E_id: gen double `ll1'=`p1'*exp(sum(log(`l1')))
			by $S_E_id: gen double `l2'=`lU2'+`lE2' +`lU2LF'+`lE2LF'
			by $S_E_id: gen double `ll2'=`p2'*exp(sum(log(`l2')))




			by $S_E_id: gen byte `last' = (_n==_N) 
	
			* individual contribution to the log likelihood 
			gen double `lnfie' = cond( !`last'  , 0, $weight * ln(( `ll1') + ( `ll2')   ))
				
					
			mlsum `lnf' = `lnfie'

	}
end

	
	
	

set type double		


tempvar mysamp 
		tempname b b01 b02 b1 b2 b03 b04 b3 b4 lnf V  a1 a2 a3  a4 b4 masshE masslE masshSE masslSE ///
		       logitp2 logitp3 logitp4  masslSELF masslELF  masshSELF masshELF

  *now put all the variables togetherlocal
 *   se_netincdiff se_convexity se_female_1 se_ageb se_edu_2-se_edu_3  se_yr_2-se_yr_4 se_reg_2-se_reg_5
 *    lnt_se t_se t2_se 
 *    e_netincdiff e_convexity e_female_1 e_ageb e_edu_2-e_edu_3  e_yr_2-e_yr_4 e_reg_2-e_reg_5
 *  lnt_e  t_e t2_e
 *  self_netincdiff self_convexity  self_female_1 self_ageb self_edu_2-self_edu_3 self_yr_2-self_yr_4 self_reg_2-self_reg_5
 * t_self t2_self lnt_self
 * elf_netincdiff elf_convexity elf_female_1 elf_ageb elf_edu_2-elf_edu_3  elf_yr_2-elf_yr_4 elf_reg_2-elf_reg_5
 * t_elf t2_elf lnt_elf
global varse " se_netincdiff se_convexity se_female_1 se_ageb se_edu_2-se_edu_3  se_yr_2-se_yr_4 se_reg_2-se_reg_5"
global vare  " e_netincdiff e_convexity e_female_1 e_ageb e_edu_2-e_edu_3  e_yr_2-e_yr_4 e_reg_2-e_reg_5"
global varself " self_netincdiff self_convexity  self_female_1 self_ageb self_edu_2-self_edu_3 self_yr_2-self_yr_4 self_reg_2-self_reg_5"
global varelf  "elf_netincdiff elf_convexity elf_female_1 elf_ageb elf_edu_2-elf_edu_3  elf_yr_2-elf_yr_4 elf_reg_2-elf_reg_5"
global dumse "  lnt_se"
global dume "    lnt_e"
global dumself " lnt_self "
global dumelf " lnt_elf"

 
set more off	
// initial values for ll function taken from clogclog		
glm durvar  $varse  $dumse if sedum==1 & lcens_spell ==0 ,  f(b) l(cloglog) 
matrix `b01' = e(b)
matrix coleq `b01' = hazardE
	
glm durvar   $vare  $dume  if sedum==0 & lcens_spell ==0   ,  f(b) l(cloglog) 	
matrix `b02' = e(b)
matrix coleq `b02' = hazardU

glm durvar  $varself  $dumself  if sedum==1 & lcens_spell ==1 ,  f(b) l(cloglog) 
matrix `b03' = e(b)

matrix coleq `b03' = hazardELF
	
glm durvar   $varelf  $dumelf  if sedum==0 & lcens_spell ==1 ,  f(b) l(cloglog) 	
matrix `b04' = e(b)

matrix coleq `b04' = hazardULF

local LL1 = e(ll)
// create matrix with initial betas from cloglog and initial mass point location and probability (given)

matrix `a1' = (-0.1 )
matrix colnames `a1' = a1:_cons

matrix `a2' = (-1.3)
matrix colnames `a2' = a2:_cons

matrix `a3' = (-3.2)
matrix colnames `a3' = a3:_cons

matrix `b1' = (-0.2 )
matrix colnames `b1' = b1:_cons

matrix `b2' = (-1.2 )
matrix colnames `b2' = b2:_cons
matrix `b3' = (-1.5 )
matrix colnames `b3' = b3:_cons


matrix `logitp2' = (0.5)
matrix colnames `logitp2' = logitp2:_cons
matrix `logitp3' = (-2.1)
matrix colnames `logitp3' = logitp3:_cons
matrix `logitp4' = (-1.1)
matrix colnames `logitp4' = logitp4:_cons


matrix `b1' = `b01',`b02',`b03',`b04',`a1',`a2',`b1',`b2',`logitp2'


matrix list `b1'

// define macro that identify id variable
global S_E_id "id"
// define macro that identify se or E spell
global SE_E "sedum"
// define macro that identy the spell
global spell "a_spell"
// define macro that identify event
global event "durvar" 
 // define macro that identify left_censored spells
global LF "lcens_spell" 
 
// define weight
global weight "pweight"
 
 
sort id a_spell year
 

//ml model, d0, no gradient or Hessian provided
ml model d0 prog2_lnf (hazardE: durvar = $varse  $dumse ) ///
                       (hazardU: durvar =  $vare  $dume  ) ///
		       (hazardELF: durvar = $varself $dumself   ) ///
		       (hazardULF: durvar = $varelf $dumelf  ) ///
					    (a1: ) (a2: ) (b1: ) (b2: )    ///
					    (logitp2: )   ///
                                , waldtest(0) init(`b1')   
/*
ml init hazardE:se_netincdiff    		=  -.4311149 
ml init hazardE:se_convexity    		=   .0589467 
ml init hazardE:se_female_1    			=  -.0197083 
ml init hazardE:se_ageb    				=  -.0122189 
ml init hazardE:se_edu_2    			=   .0085957 
ml init hazardE:se_edu_3    			=   .2173636 
ml init hazardE:se_yr_2    				=  -.0847899 
ml init hazardE:se_yr_3    				=   .0280419 
ml init hazardE:se_yr_4    				=  -.2040635 
ml init hazardE:se_reg_2    			=   .0395028 
ml init hazardE:se_reg_3    			=   .0180292 
ml init hazardE:se_reg_4    			=   .0412137 
ml init hazardE:se_reg_5    			=   .0964456 
ml init hazardE:lnt_se    				=  -.5173293 
ml init hazardE:_cons    				=  -1.150816 
ml init hazardU:e_netincdiff    		=   1.512442 
ml init hazardU:e_convexity    			=  -.1542051 
ml init hazardU:e_female_1    			=   .6027568 
ml init hazardU:e_ageb    				=   .0295773 
ml init hazardU:e_edu_2    				=   .1027266 
ml init hazardU:e_edu_3    				=   .1034038 
ml init hazardU:e_yr_2    				=   .1936451 
ml init hazardU:e_yr_3    				=   .2660251 
ml init hazardU:e_yr_4    				=  -.1771764 
ml init hazardU:e_reg_2    				=  -.2091639 
ml init hazardU:e_reg_3    				=   -.318311 
ml init hazardU:e_reg_4    				=  -.4118783 
ml init hazardU:e_reg_5    				=  -.4677638 
ml init hazardU:lnt_e    				=  -.4903883 
ml init hazardU:_cons    				=  -3.250502 
ml init hazardELF:self_netincdiff    	=  -.7752243 
ml init hazardELF:self_convexity    	=  -.0054885 
ml init hazardELF:self_female_1    		=   .1707384 
ml init hazardELF:self_ageb    			=   -.034369 
ml init hazardELF:self_edu_2    		=  -.1368612 
ml init hazardELF:self_edu_3    		=   .0442427 
ml init hazardELF:self_yr_2    			=  -.1578914 
ml init hazardELF:self_yr_3    			=  -.2990763 
ml init hazardELF:self_yr_4    			=  -.7121819 
ml init hazardELF:self_reg_2    		=    .038384 
ml init hazardELF:self_reg_3    		=  -.0402288 
ml init hazardELF:self_reg_4    		=   .1816011 
ml init hazardELF:self_reg_5    		=   .1395922 
ml init hazardELF:lnt_self    			=  -.0165969 
ml init hazardELF:_cons    				=  -.8937651 
ml init hazardULF:elf_netincdiff    	=   1.650138 
ml init hazardULF:elf_convexity    		=  -.1679822 
ml init hazardULF:elf_female_1    		=   .8080704 
ml init hazardULF:elf_ageb    			=  -.0467125 
ml init hazardULF:elf_edu_2    			=  -.0180253 
ml init hazardULF:elf_edu_3    			=   .1218535 
ml init hazardULF:elf_yr_2    			=   .1180691 
ml init hazardULF:elf_yr_3    			=   .0018832 
ml init hazardULF:elf_yr_4    			=  -.7976231 
ml init hazardULF:elf_reg_2    			=  -.2944714 
ml init hazardULF:elf_reg_3    			=  -.4856471 
ml init hazardULF:elf_reg_4    			=  -.5629266 
ml init hazardULF:elf_reg_5    			=  -.5293553 
ml init hazardULF:lnt_elf    			=  -.2430833 
ml init hazardULF:_cons    				=  -2.006452 
ml init a1:_cons    					=  -.532556 
ml init a2:_cons    					=  -1.330282  
ml init b1:_cons    					=  -3.098831 
ml init b2:_cons    					=  -1.836795 
ml init logitp2:_cons    				=  -1.415653*/
					
ml maximize, trace  search(off) difficult gradient

eststo m
eststo r1lr: nlcom (p1: 1/(1+exp([logitp2]_cons))) ///
                   (p2:exp([logitp2]_cons)/(1+exp([logitp2]_cons))), post

log close
